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Abstract 



The stochastic dynamics of inclusions in a randomly fluctuating biomembrane is sim- 
ulated. These inclusions can represent the embedded proteins and the external particles 
arriving at a cell membrane. The energetics of the biomembrane is modelled via the 
Canham-Helfrich Hamiltonian. The contributions of both the bending elastic-curvature 
energy and the surface tension of the biomembrane are taken into account. The biomem- 
brane is treated as a two-dimensional sheet whose height variations from a reference frame 
is treated as a stochastic Wiener process. The lateral diffusion parameter associated with 
this Wiener process coupled with the longitudinal diffusion parameter obtained from the 
standard Einsteinian diffusion theory completely determine the stochastic motion of the 
inclusions. It is shown that the presence of surface tension significantly affects the over- 
all dynamics of the inclusions, particularly the rate of capture of the external inclusions, 
such as drug particles, at the site of the embedded inclusions, such as the embedded 
proteins. 

PACS: 87.20; 34.20; 87.22BT 
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1 Introduction 



Amphiphilic molecules, such as lipids and proteins, are composed of hydrophilic heads and 
hydrophobic hydrocarbon-based chains. They can self-assemble themselves into a variety of 
exotic structures when placed in an aqueous environment [1]. One such structure is the bio- 
logical membrane composed of a bilayer of phospholipid molecules. Biomembranes have very 
interesting spatial dimensions. Their thicknesses can vary up to few nano-mctcrs, but their lin- 
ear sizes can take up to tens of micro-meters. They can, therefore, be regarded as highly flexible 
fluid-like two-dimensional (2D) sheets embedded in three-dimensional embedding space. 

Biomembranes form large encapsulating bags, called vesicles, because open sheet-like con- 
figurations would involve a large energy along the hydrophobic edges [2]. They play a very 
significant role in many of the life's processes, since not only they act as barriers to main- 
tain the structural integrity of cells, but also provide the functional environment for many of 
the embedded proteins that penetrate the biomembrane thicknesses and act as ion transport 
channels to the interior of the cells [3] . 

Biomembranes can be studied from three perspectives [2] , namely (i) their molecular archi- 
tecture, focusing on their material properties, (ii) as statistical-mechanical systems displaying 
a very rich and exotic variety of configurations as 2D surfaces, and (iii) their functioning in 
biological systems. 

In our previous work [4], as well as in this paper, we have been interested in biomembranes 
from the second perspective (ii) , and in particular we have been concerned with computational 
modelling of the diffusion dynamics of objects (inclusions) existing within and on the surface 
of stochastically fluctuating biomembranes. The stochastic fluctuations of the biomembranes 
can originate from thermal fluctuations. They promote shape fluctuations and shape transfor- 
mations in biomembranes. One such shape transformation, observed under video microscopy, 
is the so-called budding transition in which a change of temperature from T — 27.2 C to T = 41 
C resulted in a change of shape of a spherical vesicle into a prolate ellipsoid, and the eventual 
fission of the biomembrane [2]. 

Two types of inclusion can be considered in a biomembrane. These are the internal inclu- 
sions, or the biomembrane's own inclusions, which we refer to as the M-type, and the external 
inclusions, which we refer to as the S-type. The S-type inclusions, such as drug particles, arrive 
at the biomembrane from outside. The M-type inclusions are the embedded objects, such as 
proteins, or ion channels, that are an integral part of a biomembrane. Although embedded 
in the amphiphilic environment, these inclusions are mobile and can freely diffuse across the 
biomembrane. 

The M-type inclusions penetrate the thickness of the biomembrane, and hence disturb its 
geometry locally by forcing the amphiphilic bilayer to adjust its thickness locally so as to match 
its hydrophobic regions with those of the inclusions [5, 6], as shown in Fig.l. This change 
in the local geometry generates perturbations in the fluid-like biomembrane structure and 
these perturbations, in turn, lead to both short- and long-range biomembrane-induced forces 
between the inclusions. Such geometro-dynammic forces affect the behaviour of inclusions in a 
biomembrane in addition to any other forces, such as direct molecular forces, that may operate 
between these inclusions, or between them and the embedding biomembrane. If these induced 
perturbations are of long wavelength, then they generate long-range forces, otherwise, they 
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lead to short-range forces that act in the immediate vicinity of the inclusions. 

The deformation in local geometry is one of the three deformation modes of a biomembrane. 
Since a biomembrane is surrounded by solutions, such as water, the surface tension existing at 
the biomcmbrane-water interface changes the overall surface area of the biomembrane. This 
change constitutes the second mode of deformation. The third mode is associated with the 
bending property of a biomembrane, which is associated with its elastic-curvature property. It 
is this property that distinguishes a fluid-like biomembrane sheet from a simple fluid sheet in 
which surface tension alone dominates the surface dynamics. [2, 7]. 

In the studies concerned with the inclusion-induced local deformations, two approaches have 
been adopted to model these deformations. In one model [8, 9, 10, 11, 12] the biomembrane 
energy is taken to be the sum of the energy due to molecular compression/expansion brought 
about by the change in local geometry and the energy due to the change in the overall surface 
area, while in another model [5, 13] the contribution of the bending energy is also included. 
The first model shows that the inclusion- induced deformations force a change of thickness in 
local geometry which decays exponentially from the inclusion-imposed value to the equilibrium 
thickness of the biomembrane. This is shown schematically in Fig.l for two rod- like inclusions. 
The second model shows that the contribution of bending energy could significantly affect the 
deformation profile at the inclusion boundary, as well as the biomembrane-induced interactions. 

The S-type inclusions arrive from the outside of the biomembrane and reside on its surface. 
Fig. 2 represents a schematic representation of two such inclusions lying on the surface of a 
biomembrane sheet. 

In our previous work [4] , we performed a computational simulation of the diffusion dynamics 
of both the external and internal inclusions in a stochastically fluctuating biomembrane in 
which the biomembrane energy was solely described in terms of the bending elastic-curvature 
energy. Furthermore, in that work we considered only the contribution of the lateral component 
of the diffusion coefficient of the inclusions, brought about by the lateral fiuctuations in the 
height of the biomembrane, to the motion of the inclusions. In the present paper, we have 
extended our previous modelling on two levels. Firstly, have also included the contribution 
of biomembrane's surface tension to the stochastic dynamics of inclusions, and secondly, the 
contribution of the longitudinal component of the diffusion coefficient has also been included. 
The combination of these two components now forms the diffusion coefficient governing the 
motion of the inclusions. Various forms of this combination has been tested. It is seen that 
the dynamics of inclusions is significantly affected by the presence of surface tension. 

2 Energetics of the biomembrane in presence of surface 
tension 

Biomembranes are characterized by their bending rigidity [14, 15, 16] and they strongly un- 
dulate under thermal fluctuations. Consequently, in studying biomembranes, surface tension 
effects are usually neglected and only the bending rigidity is considered as the source control- 
ling their shape and fiuctuations. However, surface tension can contribute significantly to the 
energy of a biomembrane if geometrical constraints [17] or external perturbations, due to the 
attached proteins in biological systems, are introduced [18]. For instance, it has been shown 
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that [19] the area dependence of surface tension can lead to a stable hole of finite size in the 
biomembrane. Furthermore, interesting phenomena have been observed when biomembranes 
are perturbed by the action of microscopes [20], or in recent years by laser tweezers [21]. The 
typical effect of such tweezers on an initially flaccid biomembrane has been reported [21] where 
the application of the laser tweezers on a vesicle over few minutes caused the shape fluctuations 
of the biomembrane to disappear. With this disappearance, the biomembrane became quite 
taut, indicating the presence of a non-negligible surface tension, a > 10~^erg/cm^. This tension 
caused the expulsion of a vesicle from inside of another vesicle. Such an spontaneous expulsion 
is quite surprising [19], since it involves the opening of a large hole, of energy much larger than 
the thermal energy. In comparison, using mechanical manipulation, rather than tweezers, of 
very large vesicles, tensions as low as, cr > 10~^erg/cm^ [20] have also been recorded. It is, 
therefore, quite reasonable to assume that surface tension in the above range plays a significant 
role in the dynamic behaviour of a biomembrane. These experiments, therefore, show that a 
realistic study of the energetics of biomembranes should include not only the bending rigidity, 
but also the effect of the surface tension. 

The free elastic energy of a symmetric, nearly fiat, biomembrane sheet is described by the 
Canham-Helfrich Hamiltonian [22, 23] 

J d^a^ {(To + 2kH^ + kK] , (1) 

where 

K = det{Kij) ^ 



R1R2 

1„ . 1/1 1 



are the Gaussian and mean curvatures of the sheet respectively, Ri and R2 are the principal 
radii of the sheet, (Jq is the surface tension, k, is the bending rigidity, R is the Gaussian rigidity, 
g is the determinant of the metric tensor and a = (ai, £12) is the local coordinate defined on the 
sheet. The last term on the RHS of (1) is, by Gauss-Bonnet theorem, a topological invariant 
and does not affect the dynamics of the biomembrane if there were no changes in its topology. 
We will, therefore, ignore this term as we assume that no such changes occur. We thus have 

j Sa^ {ao + 2kH''] . (2) 

To facilitate the study of a nearly flat biomembrane, whose free energy is given by (2), it is 
convenient to consider it to be parallel to the {xi — X2) plane, regarded as the reference plane. 
The position of a point on the biomembrane can then be described by a single-valued height 
function, h, representing the position of a point on a fluctuating, nearly flat, sheet relative to 
the reference plane. In this, the so-called Monge representation [2], equation (2) is written as 



(3) 



In this form, the Hamiltonian is expressed solely in terms of the height function, h, and its 
derivatives. We remark that in our previous work [4] the contribution of the second term on 
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the RHS of (3) was not included as we did not consider the contribution of the surface tension 
ctq. The Hamiltonian in (3) describes the energetics of the biomembrane from which we obtain 
the stochastic lateral motion of both the M-type and the S-type inclusions. 

2.1 Computation of the lateral diffusion coefficient of inclusions 

Let us first consider the lateral stochastic motion of the inclusions, i.e. that component of 
the motion which is associated with the changes in h{xi,X2). One way to introduce stochastic 
behaviour into the dynamics of the biomembrane is to treat h, obtained from (3), as a stochas- 
tically fluctuating Wiener process. Evidently, this stochastic behaviour is communicated to 
the inclusions residing inside, as well as on the surface, of the biomembrane. 

To derive the lateral diffusion coefficient, D±, associated with the random changes in 
the height function, we first need to obtain the the height-height correlation function of the 
biomembrane in real space. This is obtained from (3) by first going over to the Fourier repre- 
sentation of h, i.e. 

Mx) = / (0^Mq)e''»-'' , (4) 

where x = {xi,X2) is the vector position of the point on the biomembrane. Consequently, (3) 
now reads 

= ^^"^^ ^ """^'^ h{q)h*{q) , (5) 

where * denotes the complex conjugate. The static height- height correlation function can be 
calculated, with the result 

(Mq; 0)/i*(q'; 0)) = ( Jf ) (27r)^5(q - q'), (6) 

where the averaging is done with respect to the Boltzmann weight factor ex^{—l-i/k^T)^ with 
being the Boltzmann constant. The corresponding dynamic correlation function can be 
written as 

Mq;i)=Mq;0)e-^('')*, 

leading to 

(Mq;t)/.-(q';0)> = e-*..(-^ 

where the damping factor, 7(g), reflecting the long-range character of the hydrodynamic damp- 
ing, is given by [2] 

7(?) = J- ('^?' + <7m) ■ (8) 



(27r)^5(q-q'), 



(7) 
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Here, rj denotes the coefficient of viscosity of the fluid surrounding the biomembrane. We note 
that in the absence of ao, (8) reduces to the damping factor given in equation (13) of our 
previous work [4]. In real space, the Fourier transform of (7) is given by 



(/i(x;t)/i(x;0)) 



// 



(27r)2 (27r)= 



-(/i(q;t)r(q';0))e-''-('i-'i'). 



(9) 



or substituting from (7) and (8) 




e 



+ (Jo? 



(10) 



where L is the linear size of the biomembrane and a is a molecular cut-off, of the order of 
nanometers, and /3 = The Fourier transform expression in (10) is quite complicated to 
perform and it was computed in the Mathematica software [24] by expanding the integrand 
and computing it term by term. The final analytical result, consisting of many terms, need not 
be quoted here and can be obtained from the authors upon request. In our previous work [4], 
the Fourier transform to real space was performed under simplifying assumptions regarding 
the transform of the damping factor term 7. In this paper, however, no assumptions were 
applied and the transform was fully computed. 

Now, treating h as a stochastic Wiener process implies that the mean and variance are 
given by 



in which D±, the diffusion constant associated with the random fluctuations in the height 
function at the local position {xi, X2), represents a measure of these fluctuations. Such random 
fluctuations can cause a roughening of the biomembrane surface as has been observed in NMR 

experiments [25]. 

Consequently, since the biomembrane's stochastic behaviour is transmitted to the inclu- 
sions, we can make the assumption that the centre of mass of an inclusion coinciding with the 
point {xi,X2) on the biomembrane would also experience the same lateral fluctuations, whose 
measure is the same diffusion coefficient, D±. 

Therefore, the expression for Dj_ for the inclusions in presence of surface tension can be 
obtained by comparing (10) and (12), leading to 



2.2 Computation of the longitudinal diffusion coefficient of inclu- 



In addition to the lateral component of the diffusion coefficient, computed above, we also 
require the longitudinal component, D\\, for both the S-type and the M-type inclusions. This 



{h{xi,X2;t)) = 0, 
{h(xi,X2;t)h(xi,X2;0)) = ^2D^t, 



(11) 
(12) 




(13) 



sions 
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coefficient is calculated via Einstein's general model of Brownian dynamics [26] according to 
which 

1 ^ 

Dii - Tiit) - rm \') , (14) 

where t is the delay (correlation) time, (■ ■ ■) refers to averaging over time, s is the dimensionality 
of diffusion space (2 in the present case) and N is the number of inclusions. In the actual 
computation of D\\, the time average is replaced with a summation over time origins from 
which the delay time is measured, so that (14) is written as [26] 



^ X N 

D\\ = E E [I nitk + t)- r,{tk) \'] , (15) 



k=l 1=1 

where X is total number of available time origins. 



2.3 Computation of the stochastic trajectories 

In our present numerical simulations, just as in the previous study [4] , the equations of motion 
of both the M-type and the S-type inclusions are described via the Ito stochastic calculus [28] 
whose stochastic differential equation is given by 

dr{t) = A[r(i), t] + D^^dWit). (16) 

This equation describes the random trajectories, r(t), of the centres of mass of the inclusions 
in terms of a dynamical variable of the inclusions, A[r(t),t], which is referred to as the drift 
velocity, and a term, dW{t), which is modelled by a Wiener process with the mean and variance 
given by 

{dW{t)) = 0, 
{dWi{t)dWj{t)) = 26ijdt. (17) 

In (16), Dt represents the total diffusion coefficient of the inclusions. In this paper, following 
the prescription given by [27] for the Brownian motion of rod-like objects, we represent this 
coefficient by a linear combination of the lateral and longitudinal coefficients and, to test the 
various possibilities, consider the following different combinations 

3 3" 

= Id, + Id,. (18) 

In the vicinity of an M-type inclusion, care should be exercised in using Df. Near an 
M-type inclusion, localized exponentially-decaying deformation would occur in the membrane 
geometry [4, 5], as depicted in Fig.l. This implies that near such an inclusion we can scale 
the diffusion coefficient by an exponentially decaying function, i.e. if Tq is the radius of an 
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M-type inclusion, then within a circular region R + tq, the modified diffusion coefficient could 
be modelled as 

D^t = Ae^ , (19) 

implying that when an S-type inclusion enters a region centered around an M-type inclusion, 
within the radius R+tq, the diffusion coefficient goes over to D^t and, to a first approximation, 
this is how an M-type inclusion interacts with an S-type one. 

The Ito equation (16) predicts the increment in position for a meso-scale time interval, 
dt, as a combination of deterministic and stochastic diffusive parts represented by the terms 
A[r(t)] and Dl^'^dW{t) respectively. This equation resembles the 'position' Langevin equation 
describing the Brownian motion of a particle [29]. The position Langevin equation corresponds 
to the long-time (diffusive time) configurational dynamics of a stochastic particle in which 
its momentum coordinates are in thermal equilibrium and hence can be removed from the 
equations of motion. Since we are interested in diffusive time scales as well, we can re-write 
(16) as 

dr{t)^^F{t) + Dy'dW{t), (20) 

where F(t) is the instantaneous systematic (Newtonian) force experienced by the centre of 
mass of the ith inclusion. This force can obtained from the inter-inclusion potentials, V{Rij), 
according to 

F. = -i:vR.m,). (21) 

Pertinent potentials have been constructed employing a statistical mechanics based on the 
Canham-Helfrich Hamiltonian. For the M-M interaction, this potential is given by [30, 31] 

V^iRij) = -ksT^-^, (22) 

describing the biomembrane-mediated temperature-dependent long-range interactions between 
a pair of disk shape M-type inclusions that can freely tilt with respect to each other. For the 
S-S interaction, the potential is given by [32] 

VZ{R^„ e,, ej) = -ksT^^ cos'm + (23) 

where A = vrrQ is the area of an M-type inclusion of radius tq, Rij is the distance between the 
centres of mass of two inclusions i and j, and Lj are the lengths of two S-type inclusions 
making angles 6i and dj respectively with the line joining their centres of mass (Fig. 2) and T is 
the biomembrane temperature. These expressions are derived for the rod-like inclusions that 
are assumed to be much more rigid than the ambient biomembrane so that these inclusions 
cannot move coherently with the biomembrane. 
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We performed numerical simulations of the space-time motion of inclusions described by 
(20) by adopting the following iterative scheme [33] 



X{t + dt) = X{t) + -^Fx{t)dt + J2DtdtR^, 

KbJ- 



Y(t + dt) = Y(t) + ■^FY(t)dt + j2Dtdt R^, 
ksT ^ 



(24) 



where i?^ and Ry are standard random Gaussian variables chosen separately and indepen- 
dently for each inclusion according to the procedure given in [29], and Fx, Fy are the X and 
Y components of the force F. For the S-type inclusions, we treated the angles in (23) as 
independent stochastic variables described by 

e{d + dt) = e{t) + ^^^r{t)dt + ^^2Dtdte'', (25) 

where r is the torque experienced by an S-type inclusion and is given by 

V- dV'^{Rij,9i,ej) 

'^~h ' ^ ^ 

and 6'^ is the angular analogue of R% and Ry- 



3 Results and discussion 

In our simulations, we used a square membrane patch together with the following data: 



L = 


40 /^m 


K = 


10-19 J 


T) = 


10^^ J sm-3 


(To = 


5 X 10-^erg/cm^ to 5 x 10"^erg/cm^ 


T = 


300K 


m — 




Li = 


0.1 /im 


To = 


0.01 urn 



where m is the mass of an inclusion, ro is the radius of an M-type inclusion, Li is the length of 
an S-type inclusion, k is the bending rigidity and r] is the coefficient of viscosity. The range of 
values for the surface tension was chosen to be in line with the discussions presented previously 
[20, 21]. 

To implement the equations of motion, (24) and (25), we first required the computation of 
Dt- The computation of D± part of Dt was based on (13) for a correlation time of t = 10-^s. 
Table I shows the variation of the D± for both the M-type and S-type inclusions for different 
values of the surface tension, ag. It can clearly be seen that the smaller is the surface tension, 
the larger is the value of the lateral diffusion constant. 



10 



Figure 1: Two rod-like embedded (M- type) inclusions vertically inserted in a fluid biomembrane showing 
exponentially decaying thickness-matching constraints on the bilayer at the boundaries of the inclusions. Heavy 
solid lines represent amphiphilic molecules. Figure based on [6]. 




Figure 2: Two rod-like surface (S-type) inclusions placed on the surface of the biomembrane. The rods have 
lengths Li and L^, widths ei and £2 and making angles Ox and 6*2 with the line joining their centres of mass. 
Figure based on [32]. 
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Figure 3: Stochastic trajectories of a sample of S-type inclusions. Both drift and diffusive motions are clearly 
visible. 



T(K) 


o-o(J//i m^) 


D^X''Xm2s-i) 


300 
300 
300 


5x10-23 
5x10-21 
5x10-1^ 


0.25x10-^ 
0.71x10-9 
0.46x10-^° 



Table I: Computed values of the lateral diffusion constant for different values of surface tension. 

To compute D\\ on the basis of (15), we first performed a constant-temperature molecular 
dynamic (MD)-like calculation [29] in which the inclusions interacted with each other via the 
potentials given in (22) and (23). These potentials describe the M-M and S-S type interactions 
respectively. For the M-S interactions, we adopted the simple mixing rule of an arithmetic 
average of the M-M and S-S potentials. The temperature was kept constant at T = 300 K by 
applying the Nose-Hoover thermostat (heat bath) [33], and the simulation time-step was set 
at dt = 10~^s. During this simulation our aim was to obtain the positions of the inclusions as 
a function of simulation time, and these were recorded at specified time intervals, and the data 
were then used in (15) in a separate simulation code. The same correlation time, t — 10"^ s, 
that was used for the computation of D± was used for the computation of D\\. In this case, this 
correlation time corresponded to the positions of inclusions recorded after 10^ time steps. The 
value of for the M- and S-type inclusions corresponding to this delay time was obtained to 
be 



D r"' = 0.99 x 10' 



-11 



m^ s 



(28) 

Comparison of (28) and the values fisted in Table I shows that the main contribution to 
the Dt comes from the lateral component, which was also found to be comparable to the 
corresponding value of D ^ 4.4 x 10"^ m2s-^ obtained from an MD-based simulation of a 
fully hydrated phospholipid dipalmitoylphosphatidylcholine (DPPC) bilayer diffusing in the 
;2-direction [34]. The choice of correlation time, i.e. t = 10~^s, over which the values of the 
components of were computed has to be justified. We justify this choice by noticing that 
the diffusion time scale, to, for a stochastic particle of mass m is given by [27] 



mD 



(29) 



12 




ffij m ,« 



Figure 4: A series of snapshots from dynamic simulations showing the capture of the S-type inclusions (black 
spheres) at the cite of the M-typc inclusions (white spheres) for dynamic M-type inclusions. The top row, 
(a)-(c), corresponds to the value of the surface tension cr=5x 10~^'^J//um^, the middle row, (d)-(f), to the value 
(7=5 xlO~^^J//xm^, and the bottom row, (g)-(i), to the value cr=5x 10~^^J//Ltm^. Only the centres of mass of 
the inclusions are shown. 

Since the relevant time scale for Brownian particles is at least, to = 10~^ s, then in order for 
the criterion of long-time dynamics, as embodied in (20), to be applicable to our model, the 
correlation time, t, in (13) and (15) over which the diffusion constants were calculated, must 
satisfy the condition [27] 

t > to. (30) 

For the data given in (27) and the values of Df as computed from any of the relationships in 
given (18), the values of tu obtained from (29) turn out to be much smaller than t indicating 
setting the value of the correlation time at t = 10~^ was quite reasonable. 

A set of dynamic simulations, based on (24) and (25), were performed for the different 
combinations given in (18) to compute the stochastic trajectories for both the S-type and 
M-type inclusions, and the rate of capture of the S-type inclusions at the sites of the M-type 
inclusions. The temperature was set at T=300 K . The simulation time step, dt, in (24) and 
(25), was set at dt = 10^^ s, and each simulation was performed for a total time of 4000 fis, 
i.e. 4 X 10^ iterations. The total number of inclusions consisted of 13 S-type and 23 M-type. 
Table II shows the variation of the number of S-type inclusions captured at the site of M-type 
inclusions for different values of the surface tension and different combinations of Dt, with the 
values of Dj^ and D\\ given in Table I and (28). The data in this table refer to the capture of 
the S-type inclusions when both the S-type and the M-type inclusions were mobile. 



13 



Simulation data 


captured S-type 


Dt^^{2D^ + D\\) 
captured S-type 


A = i(i^± + 2L'||) 
captured S-type 




= 5 X 10-23^ T = 


300 


10 


11 


8 


ctq 


= 5 X 10-21, T = 


300 


3 


4 


3 




= 5 X 10-l^ T = 


300 


1 


1 


1 



TablcII: Variations of the mimbcr of captured S-typc inclusions with surface tension and temperature for 
different different combinations of the lateral (-D_l) and longitudinal diffusion coefficients, uo is in units 
of J//im^, and the Ds are in units of m^ s"^. Both the S-type and the M-type inclusions were mobile. 



To make a comparison with the case when the M-type inclusions were static and only the 
S-type inclusions were mobile, the variation of the number of captured S-type inclusions for 
different combinations of Df was computed, and the results for a particular value of ctq are 
listed in Table III. 



Simulation data 


captured S-tj'po 


captured S-t}'po 


A = :^P± + 2D||) 

captured S-tj'pe 


(70 = 5 X 10-21, T = 300 


3 


3 


2 



Table III: Variations of the number of adsorbed S-type inclusions for different combinations of the lateral 
{D±) and longitudinal (Dy) diffusion coefficients when only the S-type inclusions were mobile. 



Fig. 3 shows the stochastic trajectories of a set of six S-type inclusions obtained from (24) and 
(25). Both drift motion, represented by the second term on the RHS of (24), and random 
diffusive, motion coming from membrane fluctuations, superimposed on the drift motion are 
visible. 

Fig. 4 shows the snapshots of a small portion of the biomcmbranc for a simulation in which 
both the M-type inclusions and the S-type inclusions were allowed to diffuse. In this figure, 
the black spheres represent the S-type inclusions and the white spheres the M-type inclusions. 
The spheres represent the centres of mass of the inclusions. The snapshots (a-c) refer to the 
case with ctq = 5 x lO-^^J/^m^, (d-f) to the case with ao = 5 x lO-^^J/^um^, and (g-i) to the 
case with (Tq = 5 x lO^i^J/zim^. In each row of this figure, the first snapshot, c.g (a), refers to 
the initial state, the second snapshot, e.g. (b), to the state after 2000 /xs and the last snapshot, 
e.g. (c), to the state after 4000 fxs. In the initial state, e.g. (a), in this figure, the outer M-type 
inclusions were distributed regularly and the inner ones were distributed at random. For the 
S-type inclusions, these were all distributed completely at random. As can be clearly seen from 
the snapshots in Fig. 4, the smaller the value of the surface tension, i.e. the top row snapshots, 
the higher the rate of capture of the S-type inclusions by the M-type inclusions. Furthermore, 
we observed that for the higher value of the surface tension, i.e. the bottom row, there was a 
clear collective motion of the S-type inclusions, i.e. the surface tension constrained the motion 
of the individual inclusions to an organised type of motion. 

The bottom-row snapshots showed that neither type of inclusions made any significant 
diffusion during the simulation time, i.e. for this value of the surface tension, the stiffness 
of the patch practically prevented any appreciable movement of the inclusions. This implied 
that the biomembrane surface behaved in a more stiff manner than the case of low surface 
tension or the absence of surface tension as reported in [4]. The last snapshot in the first 
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row of this figure shows that the majority of the S-type inclusions were captured by the M- 
type inclusions, and that there was a total disorder in the pattern of the M- type inclusions. 
This is similar to the results in our previous simulation in the absence of surface tension [4], 
although a bigger number of inclusions were captured during that simulation. This implies 
that as the surface tension becomes smaller and smaller, the biomembrane behaves more like 
a fluid, hence allowing an easier diffusion for any type of inclusions. This can have direct 
practical implications, namely, for the passage of external particles into the interior of a cell, 
the magnitude of the surface tension can play a very crucial role. A flexible biomembrane 
allows more external particles to arrive at the site of the embedded proteins (i.e. the M-type 
inclusions), hence increasing the rate of transport of material to the interior of the cell. This 
conclusion can have useful implications for the design artificial membranes, smart drugs and 
a better drug delivery mechanism. 

To sum up, our computational simulation has been able to reveal the essential role played 
by the surface tension in the overall dynamics of the biomembrane. This is compatible with 
the experimental finding [21] in which it was found that for a surface tension value of cr > 
5 X 10~^^J//xm^ the vesicle becomes quite taut inducing an organised type of motion to the 
contents of the vesicle. In our previous simulation [4] it was shown that the rate of capture of 
the S-type inclusions was clearly affected by an increase in the temperature of the biomembrane. 
We can now add that this rate is also affected strongly by the presence of the surface tension 
in the biomembrane. 

Another finding from our simulation is that the contribution of the lateral component of 
the diffusion coefficient, D±, computed from the stochastic fluctuations of the height function 
of the biomembrane, far outweighs that of longitudinal coefficient, D\\. In view of the fact that 
the larger is the value of the surface tension, the smaller the value of the diffusion coefficient 
becomes, then if the contribution of surface tension had also been included in the computation 
of the longitudinal component of the diffusion coefficient, this would have further reduced the 
value of this component and made its relative importance, vis-a-vise D±, even less important. 
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